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Abstract 

The phase transition in a 3D array of classical anharmonic oscillators 
with harmonic nearest-neighbour coupling (discrete (j) 4 model) is studied 
by Monte Carlo (MC) simulations and by analytical methods. The model 
allows to choose a single dimensionless parameter a determining completely 
the behaviour of the system. Changing a from to +oo allows to go 
continuously from the displacive to the order-disorder limit. We calculate 
the transition temperature T c and the temperature dependence of the order 
parameter down to T = for a wide range of the parameter a. The T c from 
MC calculations shows an excellent agreement with the known asymptotic 
values for small and large a. The obtained MC results are further compared 
with predictions of the mean-field and independent-mode approximations 
as well as with predictions of our own approximation scheme. 

In this approximation, we introduce an auxiliary system, which yields 
approximately the same temperature behaviour of the order parameter, 
but allows the decoupling of the phonon modes. 

Our approximation gives the value of T c within an error of 5% and 
satisfactorily describes the temperature dependence of the order parameter 
for all values of a. 



1. INTRODUCTION 



One of the basic classification schemes for structural phase transitions consists in 
assigning it to the order- disorder or the displacive type. The displacive transition 
can be described as a freezing of a phonon mode, which shows "critical softening" 
at the phase transition point. The occurrence of a soft mode is often used as 
criterion for a displacive transition in a real systems, since the frequency of the 
phonon modes is accessible by spectroscopic experiments. 

In the order-disorder case, there are two or more locations for each atom 
in the unit cell. Occupation numbers for these locations are the same above the 
transition temperature, and differ below. Formally, as in the displacive case, the 
system can be described in "phonon" language. 

There is a simple model which shows that one can go from the order- 
disorder to the displacive type continuously 0. This model can be defined as 
a 3D cubic lattice of classical anharmonic ID oscillators with nearest-neighbour 
harmonic coupling P, ||, f|, |], |[ [7]: 



V = 4 E *l + f E < + % E(*n - *»0Mn, »'). (1) 

Z n 4 n A n ,n> 

A, B, and C are model parameters, the indices n and n' run over all oscillators, 
a(n,n') is equal to 1 for neighbouring particles and vanishes elsewhere. The 
system undergoes a phase transition from the higher symmetry to the lower sym- 
metry phase at a certain temperature T c for any A < 0, B > 0, C > 0, i.e. the 
statistical average of each coordinate x n takes a non-zero value rj =< x n > below 
T c and vanishes above. It is often convenient to express the potential (1) as 

V = J2 v ( x n) ~ Cj2 x nX n '0-(n,n'), (2) 

n n,n' 

with an "on-site" single particle potential 

A' B 

v (x) = —x 2 + —x 4 , A' = A + 12C . (3) 
2 4 

It is known that the behaviour of the system is governed by the ratio 

a = -A/C. (4) 



At small a > the system shows a displacive phase transition, while for 
large a the system behaves as the Ising model, which shows a typical order- 
disorder phase transition. The transition temperatures T c in the limit cases are 
known from Ising-model and self-consistent phonon calculations, to be respec- 
tively i, | g 1: 
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T c {a | 0) « 2.64C|A|/5, 
T c (a -f +00) » 9.12C|A|/£, (5) 

assuming here and further that the temperature is expressed in energy units (the 
Boltzman constant equal to 1). On the other hand, despite of the important role 
of the above model in the theory of the structural phase transitions 0, the actual 
dependence of T c (a) is not known. The results of previous molecular dynamics 
and Monte Carlo studies are collected in Figure 1. They obviously do not give a 
consistent quantitative picture. So far the analytical study was restricted to the 
mean-field approach. |2|, |10 . 



Recently, it was observed that the knowledge of the dependence T c (a) can be 
useful in the quantitative analysis of the properties of crystalline S , n 2 P2>S6 which 
has a ferroelectric phase transition showing simultaneously features typical for 
both the order- disorder and displacive type ]TT[ . 



The aim of this paper is to establish this dependence of T c (a) as well as 
the temperature behaviour of the order parameter. Let us stress that, similarly 
to some related papers ||, 0, || |7j , we are not interested here in details of critical 
behaviour in the very vicinity of the phase transition. Critical behaviour of this 
model is thoroughly described for example in the reference 0. 

The paper is organized as follows. Section 2 describes our MC simulations 
performed for a wide range of values of the parameter a. In section 3, we first 
compare the MC results with rather poor predictions of the standard decoupling 
schemes and suggest an improved self-consistent equation for the order parameter 
that allows to calculate both the transition temperature and the order parameter 
with a reasonable accuracy for all values of a. 



2. MONTE-CARLO SIMULATIONS 



For numerical simulation it is convenient to re-scale coordinates and energy units. 
This allows to reduce the potential energy ([]]) into the form 

^red = ~ XX + T J2 X n + \ J2( X n ~ X n >) 2 a{fl, n'), (6) 
n n n,n' 

with a single dimensionless parameter a = —A/C. Then the re-scaled order 
parameter at zero temperature is equal to 1 for any a > 0. 

The typical size of the array of atoms studied in our MC simulations is 
10 x 10 x 10 atoms, with periodic boundary conditions. We perform Monte-Carlo 
steps consecutively for each atom, and accept (or reject) them accordingly to 
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standard criteria. Additionally, we perform "magic" steps for the case of large a, 
when the sign of the coordinate of the given atom may flip. These steps allow the 
system to thermalize in the order- disorder limit as well. We calculate the square 
of the order parameter as the average 

^ = A r-i/2 <X 2_ X 2 >5 (7) 

where = N~ l l 2 J2 x n e %kn is the Fourier transform of x n , N is the total number 
of particles. 

For the case of an infinite slab, N^^X 2 is negligible, and (0) gives purely 
the square of the order parameter. For a finite slab, the term X\ allows to remove 
fluctuations from the high-temperature branch. 

The results of the calculations are presented in Figures 2-4. It is crucial 
to check the dependence of the results on the system size. Figure 2 presents the 
temperature dependence of X for sizes 15 x 15 x 15 and 5x5x5. It is clear, 
that change of the size of the slab affects practically only the fluctuation region 
near T c . The value of T c calculated from the fit of the dependence (see Figure 2) 
remains almost unchanged. This type of size dependence of the data is found for 
the whole range of a. 

Figure 3 presents data for rf{T) obtained for the potential (^) for different 
values of the parameter a. Note that the Landau theory yields a linear tempera- 
ture dependence for rj 2 (T) . 

Values of T c are extracted from the data presented in Figure 2. The plot 
for T c (a) is given in Figure 4 where a logarithmic scale for the a- axis is used. The 
monotonic dependence approaches known limit values with a good accuracy. The 
change in a is for which T c (a) varies significantly is about 2 orders of magnitude. 



3. ANALYTICAL APPROACHES 

Two standard decoupling schemes have been used in the literature to make the 
phase transition in the model tractable, usually referred to as mean field (or inde- 
pendent site) approximation and self-consistent phonon (or independent mode) 
approximation. In this section we first analyse the advantages and disadvantages 
of these standard approximations and then we propose a modified approximation 
scheme that combines advantages of both schemes. 
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3.1 Independent mode approximation 



In the independent mode approximation (IMA), the deviations from the average 
value given by the order parameter 

Vn = Xn~V (8) 

are represented by Fourier coordinates Yf, = N~ 1 ^ 2 J2 n Vn^ %kn - Interaction between 
Fourier coordinates is simplified by assuming that each Fourier coordinate is 
influenced only by the average of its interactions with the other coordinates. 
This leads to an effective harmonic approximation. The order parameter in IMA 
is defined by the equation 

Ar\ + Brf + 3Br]I(T) = 0, (9) 

where the function I(T) = N" 1 J2kYkY-k is calculated from the phonon disper- 
sion relation renormalized by the given value of the order parameter and the 
thermal fluctuations. In the vicinity of the phase transition point, I(T) can be 
evaluated by assuming a " critical" phonon dispersion (with zero frequency of the 
zone-center mode): 

where k « 2.638 . Note that the stability limit T c ima = —ACk/B of the high 
temperature phase as obtained from (|) and (|T0|), provides an exact prediction 
for T c and r\ (T) in the displacive limit. However, for the order- disorder limit, the 
IMA values differ considerably from the exact values. 



3.2 Mean field approximation 

In the mean-field approximation (MFA) for the system with a harmonic coupling, 
the behaviour of the original system is modelled by an auxiliary system in which 
all direct inter-site interactions are replaced by an effective external field E, but 
the on-site anharmonicity is kept without any approximation. Taking the on-site 
potential as given by ([3|), the ensemble averages in such an auxiliary system at 
fixed external field are given by 

, =< „ >= me) - J rf (id 

J exp| — (v(x) — Ex)/T\dx 

Since at finite temperatures the gr{E) is a monotonic function, it can be inverted 
and the self-consistent equation for the order parameter in the auxiliary system 
can be written as 

E = g T \r ] ). (12) 
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The effective field E is defined as the force on x n supplied by the interaction terms 
separated in (0), assuming that the displacement of the six nearest neighbouring 
sites (or at least their sum) is frozen at the equilibrium value 77 : 

E = YlCr]. (13) 



Self-consistent solution of equations (|12|,[T3|) defines the order parameter //mfa 



in MFA. The phase transition temperature T C MFA( a ) at which tjmfa vanishes is 
shown in Figure 1. It was previously remarked by S. Aubry [0], |2j that the relative 
over-estimation of T c by MFA is almost the same (about 30 percent) for both limit 
cases (a — > +0, a — > +oo). Comparison of T Cj mfa(o) with our MC results shows 
that the discrepancy is really systematic for all intermediate cases. Although this 
error is rather large, its systematic character strongly suggests that the physics 
of the crossover is already well taken into account by the MFA. 

Let us analyse the function gr{E) describing the auxiliary ensemble of the 
uncoupled on-site oscillators in more detail. Let us stress the following points: 



1. The variation of T c with a is within MFA entirely given by the slope of the 
function gr(v) at 7/ = 0. 

2. Unlike the on-site potential v(x), the function gr{E) at finite temperature 



is a smooth monotonic odd function (see Figure 5) at any T, a JTI| • Both 
gr(E) and its inverse g? iv) can be expanded in Taylor series: 

oo oo 

9t(E) =J2x2i-i(T)E 2i -\ gr\v) = E^-i(T)rf - 1 (14) 

i=l i=l 

3. The function g^iv) can be identified with the derivative of its free energy 
F(j], T), which can thus be written in the form 

F(^T) = F(0,T) + £%^rf (15) 

8=1 

4. Obviously, the Taylor expansion coefficients of gr(E) and g^ijf) are related 
(£1X1 = 1, ^3Xi+Xi^3 — 0, etc.) This allows to express £21-1 CO i n the limit 
case of the weak anharmonicity (B << A') by expanding gr{E) in powers 
of B. With an accuracy 0(B 2 ) we obtain 

£i(T)=A'+^f> UT) = B. (16) 

In the strongly anharmonic order-disorder limit (A' < 0,T << A' 2 /B), 
expressing gr(E) via averages < x 2 >, < x A > yields 

frCO = ^ > &CO = ^2 ' • ( 17 ) 
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Finally, let us note that in the weak anharmonicity case we can solve the 
inverse problem - express the parameters of the on-site potential via the 
first two free energy coefficients £i(T), £ 3 (T). With the same accuracy as 

(0) 

A' = 6C0 - , B = UT) ■ (18) 



3.3 Combined scheme 



We have seen that the IMA predicts well the phase transition temperature in the 
displacive limit, while MFA predicts rather well its variation with a. It would 
be desirable to have an approximate equation of state for the system (p]) that 
combines the advantages of both above discussed approaches. The key idea of 
our approach is the assumption of the existence of an effective potential (with 
temperature dependent coefficients) for which the self-consistent phonon approx- 
imation gives correctly the order parameter. In determining the coefficients of 
such an effective potential, we use the properties of the free energy F(r],T) (re- 
spectively its derivative g^iv)) °f the auxiliary system of uncoupled anharmonic 
oscillators discussed above. 

More precisely, the self-consistent equation for r](T) is constructed in three 
steps, as follows: 



1. We look for an effective on-site potential of the form 

u(x) = jx 2 + jX A (19) 

where a' and (3 are defined by the expressions that appear in the above 
discussed inverse problem (|T8|): 



«' = 6-^, = (20) 

This potential obviously coincides with v(x) in the weak anharmonic limit. 

2. We introduce a function £i >e g(T, 77), which allows to write g^ 1 formally as a 
finite polynomial: 

9t X {v)=^M t ^ V)V + UTW- (21) 



These functions £i ie s(T) are used instead of £i(T) in the definitions (p0|), so 
that we have 

*' = ZiMT,v)- 5 3( I )T , > P = &n (22) 

Note that the potential u(x) still coincides with v(x) in the weak anharmonic 
limit for small rj, since Ci,es(T,rj) goes to £i(T) for 77 — > 0. 
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3. We consider eq. (Q) and replace v(x) with coefficients A' and B by an 
expression u(x) with temperature-dependent coefficients a' and (3 defined 
in eq. fl2"2"|). Then we apply IMA to this auxiliary system. The eq. ([]) then 
becomes 



£i, eff (T,77)-12C + 



r) + t 3 (T) V 3 = 0. 



(23) 



This equation is to be solved together with formulae (pi] ) and (pTT]) , defining 
^i cfj and gx, respectively. The value of ^(T), entering these equations, is 
given by the series (|14|). 



For the calculation of the phase transition temperature only, the second step 
can be omitted. It is obvious from its construction that the suggested method 
provides the same (exact) result for the T c in the displacive limit as the usual 
IMA. In the extreme order- disorder limit, the value of T c defined by (|23|) can 
be obtained analytically using (|T7D. The resulting value of T c overestimates the 
known Ising value by less then 7%. The principal advantage of the modified 
approach is that it allows to calculate the T c (and f]{T)) with the above or better 
accuracy for all values of a, as it can be seen from the comparison with our 
MC data (Figure 6). The MC result for r] 2 (T) is satisfactorily described as well 
(Figure 7). 



4. DISCUSSION 



Let us analyze the proposed model in comparison with the standard decoupling 
schemes. The latter treat the system as a gas of elementary excitations, which 
are supposed to interact weakly. The assumption of weak interaction allows to 
replace the interaction between the elementary excitations with an interaction 
with an average field. Choice of the elementary excitations as the plane waves or 
on-site oscillators yields IMA or MFA, respectively. It is clear, however, that the 
assumed weakness of the interaction is actually not realized for the general case, 
no matter what elementary excitations we choose. 

The main advantage of our approach is that it virtually replaces the real 
strongly-correlated system ([!]) with an auxiliary one, which allows decoupling. 
It is also worth noting that the theory is carried out in terms of gr{v) > which is 
always a smooth monotonic function. Moreover, gr{v) does not change drastically 
when a is varied from to +oo - the calculation of parameters £i and £3 at T c shows 
that their dimensionless values lie within the relatively narrow ranges 9. 7... 12 and 
0...4, respectively. Therefore the replacement procedure works uniformly well for 
all values of the parameters. 

It would be interesting to investigate possible extensions to higher-order 
terms. A systematic extension of the present method should contain a larger 
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number of terms in the effective on-site potential QT9D and in the expression for 
the function g^iv) m (0); an d solve the self-consistent equation for the auxiliary 
system more accurately than . 

As a simplification we can consider a purely linear auxiliary system, i.e. 
neglect the (3 term in (|l9|) and ^(T) in fl22f) . We obtain g^iv) = £i,efi{T,r))r] 
and (^) then reduces simply to the mean-field equation of state g^iv) = 12Cr). 
Therefore, the scheme proposed here can also be considered as a generalization 
of the mean field approximation. 



Our method can be applied to more complicated models for which the 
self-consistent phonon theory is exact in the weak anharmonic limit. This is 
particularly interesting for the analysis of the DIFFOUR model fl2| in which the 
additional second neighbour harmonic coupling shows a phase transition to an 
incommensurate phase for which the MC calculations are much more difficult. 



5. CONCLUSIONS 



We have studied the crossover from a displacive to an order-disorder phase tran- 
sition in the discrete 4 model with first-neighbour coupling.. The crossover is 
governed by the single parameter a. Quantitative information about T c (a) and 
r)(T, a) in this simple model may be helpful in elucidating the behaviour of some 
real crystals with phase transitions of a mixed displacive and order-disorder type. 

In terms of the dimensionless parameter a we determined the change of 
the transition temperature by Monte-Carlo calculations. These show a crossover 
from the displacive to the order- disorder limit. 

Monte Carlo calculations have shown an excellent agreement for T c in the 
two limit cases in which exact results are known. We expect that the same 
precision is obtained for the intermediate region. Thus, the presented Monte 
Carlo results can be taken as quite reliable estimates of T c (a) with a precision of 
the order of 1% and we believe that a comparable precision was achieved for the 
temperature dependence of the order parameter (except in the critical region in 
the vicinity of the phase transition) . 

We have presented an analytical approach, which goes beyond the conven- 
tional decoupling schemes. For this, we introduce the auxiliary array of oscillators 
that (i) can be treated in the independent-mode approximation and (ii) yields 
approximately the same values of T c and the order parameter, as the real system. 
The method combines the equation of state of the self-consistent phonon theory 
with the response function of the system of uncoupled anharmonic oscillators 
used in the the mean-field theory. It can be presented as a generalization of the 
mean- field scheme. 
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The analytical results for T c agree with Monte-Carlo simulations with about 
5% accuracy. Further improvement could possibly come from higher order terms 
in the expansion we have used. The formalism can be used to study incommen- 
surate phase transitions as well. 
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Figure captions 



Fig.l The critical temperature vs. In a. The curve shows the mean-field 
result; thin horizontal lines show the asymptotic values of T c . Results of previous 
Molecular-Dynamics |5|, |3[ and Monte-Carlo [|J calculations are plotted with filled 
and open circles, respectively. 

Fig. 2 The role of the finite size of the system studied numerically. Numerical 
data for the square of the order parameter plotted as a function of T. Filled 
circles: 5x5x5 oscillators; open circles: 15 x 15 x 15 oscillators. Both results are 
obtained for a = 5 in eq. (0) by averaging over 3000 realizations at each point. 
The solid line shows the interpolation used to obtain T c . 

Fig. 3 The temperature dependence of the square of the order parameter for 
values of a varying from 0.98 to 4000. There is a factor \[2 between the a values 
for the neighbouring curves. The data are obtained using eq. (|7|) by averaging 
over 1000 realizations at each point; the relative amount of "magic" steps is 0.02. 

Fig. 4 Numerical results for the critical temperature T c vs. In a. The values 
of T c are extracted from the data presented in Fig.3. Thin horizontal lines show 
the asymptotic values of T c . 

Fig. 5 Typical dependence of g^iv) & t small a (solid line) and large a 
(dashed line). The inset shows the on-site potential for both cases. (Calculated 
for (^|)with T = 5 and a = 1 and 100, respectively.) 

Fig. 6 Numerical data for T c (lna) compared with results of calculations by 



equations (E5L 1 1 , 2T|) (the solid line). The mean- field approximation is given by 



the dashed line. 

Fig. 7 The temperature dependence of the square of the order parameter at 
several values of a(= 0.98, 3.9, 15.6, 62.5, 250, 1000, 4000). T c grows with increas- 
ing a: numerical data (points) and calculation from (p3|,p4],|2T|) (lines). 
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